derivatives_dir <- "/Volumes/django/EP_EEGfMRI_deconvolution/derivatives"
output_dir <- "/Volumes/django/EP_EEGfMRI_deconvolution/derivatives/results"
# The neuRosim::canonicalHRF() implements the standard double‑gamma haemodynamic response function as originally characterized in event‑related fMRI by Friston et al. (1998) and refined by Glover (1999)
# Parameter Description Default
# a1 Delay of the main (positive) response 6s
# a2 Delay of the undershoot 12s
# b1 Dispersion (width) of the main response 0.9
# b2 Dispersion of the undershoot 0.9
# c Scale (height) of the undershoot relative to the main response 0.35
hrfs <- read_csv ('./hrflib18.csv' , show_col_types = FALSE )
hrf_pcs <- read_csv ('./HRFbasis_500Hz.csv' , show_col_types = FALSE ) %>%
rename (PC1 = pc1, PC2 = pc2, PC3 = pc3) %>%
arrange (time) %>%
pivot_longer (
cols = - time,
names_to = "pc" ,
values_to = "value"
)
# default
hrf_canonical <- canonicalHRF (hrfs$ hrf_time)
long_hrf <- hrfs %>%
arrange (hrf_time) %>%
pivot_longer (
cols = - hrf_time,
names_to = "hrf" ,
values_to = "value"
) %>%
group_by (hrf) %>%
mutate (
corr = cor (value, hrf_canonical),
rmse = sqrt (mean ((value - hrf_canonical)^ 2 ))
) %>%
ungroup () %>%
mutate (
hrf_group = case_when (
hrf == "hrf1" ~ "group2" ,
hrf == "hrf2" ~ "group2" ,
hrf == "hrf3" ~ "group1" ,
hrf == "hrf4" ~ "group2" ,
hrf == "hrf5" ~ "group1" ,
hrf == "hrf6" ~ "group2" ,
hrf == "hrf7" ~ "group1" ,
hrf == "hrf8" ~ "group1" ,
hrf == "hrf9" ~ "group1" ,
hrf == "hrf10" ~ "group3" ,
hrf == "hrf11" ~ "group3" ,
hrf == "hrf12" ~ "group3" ,
hrf == "hrf13" ~ "group4" ,
hrf == "hrf14" ~ "group3" ,
hrf == "hrf15" ~ "group4" ,
hrf == "hrf16" ~ "group4" ,
hrf == "hrf17" ~ "group4" ,
hrf == "hrf18" ~ "group4" ,
) ) %>%
mutate (panel = hrf_group) %>%
group_by (hrf_group) %>%
mutate (corr_group = mean (corr),
rmse_group = mean (rmse)) %>%
ungroup ()
cor_canonical <- long_hrf %>%
distinct (hrf_group, corr_group, rmse_group) %>%
arrange (hrf_group)
ylim_hrfs <- range (long_hrf$ value) * 1.05
panel_labels <- c (
group1 = "group1: no undershoot" ,
group2 = "group2: low undershoot" ,
group3 = "group3: moderate and early undershoot" ,
group4 = "group4: large and later undershoot"
)
long_hrf %>%
ggplot (aes (x = hrf_time, y = value, group = hrf, color = rmse)) +
geom_hline (yintercept = 0 , color = "grey" , size = 0.2 ) +
geom_line (color = "white" , size = 1.6 ) +
geom_line (size = 0.8 , alpha = 1 ) +
scale_color_gradientn (
name = "RMSE to the canonical HRF" ,
colours = c (
"#fee0d2" ,
"#fcbba1" ,
"#fc9272" ,
"#ef3b2c" ,
"#cb181d" ,
"#a50f15" ,
"#67000d"
),
guide = guide_colorbar (
title.position = "top" ,
title.hjust = 0.5 ,
barwidth = 15 ,
barheight = 1
)
) +
scale_y_continuous (
breaks = 0 ,
limits = ylim_hrfs
) +
facet_wrap (
~ panel,
ncol = 2 ,
labeller = labeller (panel = panel_labels)
) +
labs (
x = "Time (s)" ,
y = NULL
) +
theme_minimal (base_size = 16 ) +
theme (
strip.text.x = element_text (hjust = 0 , size = 18 ),
axis.title.x = element_text (size = 18 ),
axis.text.x = element_text (size = 18 ),
axis.text.y = element_text (size = 14 ),
panel.grid.major.x = element_line (color = "grey80" , size = 0.2 ),
panel.grid.major.y = element_blank (),
panel.grid.minor = element_blank (),
panel.spacing = unit (1 , "lines" ),
legend.position = "bottom" ,
legend.direction = "horizontal" ,
legend.title = element_text (size = 16 ),
legend.text = element_text (size = 16 )
) -> p
p_file1 <- tempfile (tmpdir = "./tmp" , fileext = '.png' )
agg_png (p_file1, width = 10 , height = 16 / 2 , units = "in" , res = 300 )
suppressWarnings (print (p))
invisible (dev.off ())
anno_df <- hrf_pcs %>%
distinct (pc) %>%
mutate (x = 10 , y = 0 )
hrf_pcs %>%
ggplot (aes (x = time, y = value, group = pc)) +
geom_hline (yintercept = 0 , color = "grey" , size = 0.2 ) +
geom_line (size = 0.8 , alpha = 1 ) +
scale_y_continuous (breaks = 0 ) +
facet_wrap (~ pc, ncol = 1 ) +
geom_text (
data = anno_df,
aes (x = x, y = y, label = pc),
size = 6 ,
hjust = 0
) +
labs (
x = "Time (s)" ,
y = NULL
) +
theme_minimal (base_size = 16 ) +
theme (
strip.text = element_blank (),
strip.background = element_blank (),
axis.title.x = element_text (size = 20 ),
axis.text.x = element_text (size = 20 ),
axis.text.y = element_text (size = 20 ),
panel.grid.major.x = element_line (color = "grey80" , size = 0.2 ),
panel.grid.major.y = element_blank (),
panel.grid.minor = element_blank (),
panel.spacing = unit (0.2 , "lines" ),
legend.position = "none"
) -> p1
long_hrf %>%
ggplot (aes (x = hrf_time, y = value, group = hrf, color = hrf_group)) +
geom_hline (yintercept = 0 , color = "grey" , size = 0.2 ) +
geom_line (color = "white" , size = 1.6 ) +
geom_line (size = 0.8 , alpha = 0.6 ) +
scale_color_manual (values = c ("#00468BFF" , "#00468BFF" , "#00468BFF" ,"#00468BFF" )) +
scale_y_continuous (
breaks = 0 ,
limits = ylim_hrfs
) +
labs (
x = "Time (s)" ,
y = NULL
) +
theme_minimal (base_size = 16 ) +
theme (
strip.text.x = element_text (hjust = 0 , size = 20 ),
axis.title.x = element_text (size = 20 ),
axis.text.x = element_text (size = 20 ),
axis.text.y = element_text (size = 20 ),
panel.grid.major.x = element_line (color = "grey80" , size = 0.2 ),
panel.grid.major.y = element_blank (),
panel.grid.minor = element_blank (),
panel.spacing = unit (1 , "lines" ),
legend.position = "none" ,
legend.direction = "horizontal" ,
legend.title = element_blank (),
legend.text = element_text (size = 16 )
) -> p2
p <- ggarrange (p1, p2, nrow = 1 , ncol = 2 , widths = c (1 , 1 ))
p_file2 <- tempfile (tmpdir = "./tmp" , fileext = '.png' )
agg_png (p_file2, width = 10 , height = 4.2 , units = "in" , res = 300 )
suppressWarnings (print (p))
invisible (dev.off ())
# cavity_tibble contains all voxel xyz inside the cavity per subj per surgery
cavity_tibble <- read_csv ('./cavity_tibble.csv' , show_col_types = FALSE ) %>%
rename (postop_date = surgery) %>%
mutate (postop_date = as.Date (postop_date, format = "%d_%m_%Y" )) %>%
group_by (subj, voxel) %>%
slice_min (postop_date, n = 1 , with_ties = FALSE ) %>%
ungroup () %>%
mutate (postop_date = format (postop_date, "%d_%m_%Y" ))
# hrflib contains all HRF information per voxel
hrflib <- read_csv ("./table_hrflib18_withsign.csv" , show_col_types = FALSE ) %>%
mutate (type_dur = ifelse (str_detect (type_name, "-dur" ), "long-lasting" , "spike" )) %>%
mutate (subj = str_replace (subj, "sub_" , "sub-" ))
fs_label <- read_csv ("./afni_fs_aparc+aseg_2009_goodnames.csv" , show_col_types = FALSE ) %>%
rename (fs_label = AFNI_label) %>%
mutate (lobe = str_c (hemisphere, " " , lobe))
# ep_info contains clincal information with surgery date
ep_info <- read_csv ("./ep_info.csv" , show_col_types = FALSE ) %>%
mutate (age_scan = round (age - as.numeric (difftime (surgery_date, scan_date, units = "days" )) / 365 ),
duration_of_epilepsy_scan = round (duration_of_epilepsy - as.numeric (difftime (surgery_date, scan_date, units = "days" )) / 365 ),
surgery_date = format (surgery_date, "%d_%m_%Y" ),
postop_date = format (postop_date, "%d_%m_%Y" )) %>%
rename (subj = participant_id,
surgery = surgery_date) %>%
group_by (subj) %>%
mutate (age_scan = round (mean (age_scan)),
duration_of_epilepsy_scan = round (mean (duration_of_epilepsy_scan))) %>%
ungroup () %>%
select (subj, gender, age_scan, duration_of_epilepsy_scan, pathology, surgery, postop_date, MRIlesion)
# this step match the clinical info particular to patient and surgery (by date)
cavity_tibble <- cavity_tibble %>%
left_join (ep_info, by = c ("subj" , "postop_date" )) %>%
select (- c ("postop_date" ))
# add voxel label
hrflib <- hrflib %>%
mutate (
hrf_group = case_when (
hrf_id == "1" ~ "2" ,
hrf_id == "2" ~ "2" ,
hrf_id == "3" ~ "1" ,
hrf_id == "4" ~ "2" ,
hrf_id == "5" ~ "1" ,
hrf_id == "6" ~ "2" ,
hrf_id == "7" ~ "1" ,
hrf_id == "8" ~ "1" ,
hrf_id == "9" ~ "1" ,
hrf_id == "10" ~ "3" ,
hrf_id == "11" ~ "3" ,
hrf_id == "12" ~ "3" ,
hrf_id == "13" ~ "4" ,
hrf_id == "14" ~ "3" ,
hrf_id == "15" ~ "4" ,
hrf_id == "16" ~ "4" ,
hrf_id == "17" ~ "4" ,
hrf_id == "18" ~ "4" ,
) ) %>%
left_join (fs_label %>% select (fs_label, labelname, lobe, tissue), by = "fs_label" ) %>%
left_join (cavity_tibble %>% distinct (subj, gender, age_scan), by = "subj" )
# add pathology
hrflib_clean <- hrflib %>%
left_join (cavity_tibble %>% distinct (subj, voxel, engel, pathology, surgery, MRIlesion), by = c ("subj" , "voxel" )) %>%
mutate (pathology_complex = pathology) %>%
mutate (pathology = case_when (
pathology_complex == "Cortical malformation" ~ "cortical malformation" ,
pathology_complex == "non lesional" ~ "normal" ,
pathology_complex == "HS" ~ "HS|dual" ,
pathology_complex == "dual pathology" ~ "HS|dual" ,
pathology_complex == "Gliosis" ~ "gliosis" ,
pathology_complex == "encephalitis" ~ "enceph|tumor|vascular" ,
pathology_complex == "Vascular malformation" ~ "enceph|tumor|vascular" ,
pathology_complex == "complex pathology" ~ "complex" ,
pathology_complex == "Tumor" ~ "enceph|tumor|vascular" ,
TRUE ~ NA )) %>%
rename (age = age_scan) %>%
filter (abs (r) >= 0.7 ,
! is.na (lobe),
tissue == "tiss__gm" ) %>%
filter (! fs_label %in% c ("6" , "13" , "26" )) %>%
filter (
abs (posPeakAmp) < 2.5 ,
posPeakTime > 3 & posPeakTime < 13 ,
abs (negPeakAmp) < 2.5 ,
)
df_all <- hrflib_clean %>%
group_by (subj, type, fs_label) %>%
mutate (
pathology = {
non_na_path <- pathology[! is.na (pathology)]
if (length (non_na_path) == 0 ) {
NA
} else {
most_common <- names (sort (table (non_na_path), decreasing = TRUE ))[1 ]
most_common
}
}
) %>%
ungroup () %>%
filter (is.na (pathology) | pathology == "normal" ) %>%
group_by (fs_label) %>%
mutate (voxelinparcel = n ()) %>%
ungroup () %>%
filter (voxelinparcel >= 5 ) %>%
mutate (sign = ifelse (sign == 1 , "positive" , "negative" )) %>%
mutate (gender = factor (gender),
fs_label = factor (fs_label),
labelname = factor (labelname),
subj = factor (subj),
lobe = factor (lobe),
sign = factor (sign)) %>%
mutate (hrf_id = factor (hrf_id),
hrf_group = factor (hrf_group, levels = c ("1" , "2" , "3" , "4" ))) %>%
select (subj, hrf_id, hrf_group, age, gender, type_dur, lobe, fs_label, labelname,
posPeakAmp, posPeakTime, posPeakFWHM, negPeakAmp, negPeakTime, negPeakFWHM, sign)
df <- hrflib_clean %>%
filter (! is.na (pathology)) %>%
group_by (fs_label) %>%
mutate (voxelinparcel = n ()) %>%
ungroup () %>%
filter (voxelinparcel >= 5 ) %>%
group_by (pathology_complex) %>%
mutate (voxelinparcel = n ()) %>%
ungroup () %>%
filter (voxelinparcel >= 5 ) %>%
mutate (sign = ifelse (sign == 1 , "positive" , "negative" )) %>%
mutate (pathology = ifelse (pathology == "enceph|tumor|vascular" , "encephalitis" , as.character (pathology))) %>%
mutate (gender = factor (gender),
fs_label = factor (fs_label),
labelname = factor (labelname),
subj = factor (subj),
lobe = factor (lobe),
pathology = factor (pathology),
sign = factor (sign)) %>%
mutate (hrf_id = factor (hrf_id),
hrf_group = factor (hrf_group, levels = c ("1" , "2" , "3" , "4" ))) %>%
select (subj, hrf_id, hrf_group, age, gender, type_dur, lobe, fs_label, pathology, labelname,
posPeakAmp, posPeakTime, posPeakFWHM, negPeakAmp, negPeakTime, negPeakFWHM, sign)
p_file_posPeakAmp <-
plot_histrogram (df_all, "posPeakAmp" , "BOLD % change" , title = bquote ("distribution of HRF peak amplitude" ), NULL )
p_file_posPeakTime <-
plot_histrogram (df_all, "posPeakTime" , "Time(s)" , title = bquote ("distribution of HRF time to peak" ), NULL )
p_file_posPeakFWHM <-
plot_histrogram (df_all, "posPeakFWHM" , "FWHM(s)" , title = bquote ("distribution of HRF peak FWHM" ), NULL )
p_file_negPeakAmp <-
plot_histrogram (df_all, "negPeakAmp" , "BOLD % change" , title = bquote ("distribution of HRF undershoot amplitude" ), NULL )
p_file_negPeakTime <-
plot_histrogram (df_all %>% filter (negPeakTime> 7 ), "negPeakTime" , "Time(s)" , title = bquote ("distribution of HRF time to undershoot" ), NULL )
p_file_negPeakFWHM <-
plot_histrogram (df_all, "negPeakFWHM" , "FWHM(s)" , title = bquote ("distribution of HRF undershoot FWHM" ), NULL )